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FINITE  ELEMENT  CALCULATIONS  OF  VISCOELASTIC  FLUID  FLOW  IN  A  SPINNING  AND 

NUTATING  CYLINDER 

1 .  INTRODUCTION 

In  recent  years  there  has  been  considerable  interest  in  understanding 
the  onset  of  instability  in  liquid- filled  shells  that  are  simultaneously 
spinning  and  nutating  in  flight.  Field  observations  by  D'Amico  and  Mil¬ 
ler  (1979)  indicate  that  there  is  an  instability  which  appears  in  fluids 
of  very  high  viscosity,  and  is  therefore  unrelated  to  the  instability 
associated  with  the  presence  of  inertial  waves  in  fluids  of  very  low 
viscosity  (Stewartson  1959,  Wedemeyer  1966).  Theoretical  work  by  Vaughn 
et  al  (1985)  and  Herbert  (1985),  coupled  with  the  experiments  of  Miller 
(1982)  performed  under  controlled  conditions,  confirm  the  presence  of 
this  high-viscosity  instability,  which  manifests  itself  through  a  rela¬ 
tively  large  despin  moment  (leading  to  a  loss  in  spin  rate).  It  is  gen¬ 
erally  thought  that  the  despin  moment  attains  a  maximum  at  a  fairly  high 
viscosity  value  (low  Reynolds  number)  and  then  decreases  as  the  visco¬ 
sity  decreases. 

Hitherto  all  work  on  this  problem,  both  experimental  and  theoretical, 
has  been  restricted  to  Newtonian  liquids.  There  are,  however,  practical 
reasons  that  make  It  important  to  consider  the  behavior  and  response  of 
the  spinning-nutating  system  when  the  liquid  fill  is  non -Newtonian,  and 
the  present  project  is  directed  at  comprehending  and  resolving  the 
issues  that  arise  in  that  case. 

The  primary  purpose  of  this  project  was  to  determine  the  liquid- induced 
despin  moment  fur  the  configuration  shown  in  Figure  1.  A  right  circular 
cylinder  of  length  2c  and  diameter  2a  spins  about  its  axis  with  angular 
speed  w.  The  axis  of  the  cylinder  is  inclined  to  the  vertical  at  the 


nutation  angle  9  and  rotates  about  the  vertical  with  angular  speed  0. 
The  cylinder  is  completely  filled  with  a  liquid  of  constant  density  p; 
this  liquid  is  non-Newtonian,  and  the  implications  of  this  will  be  dis¬ 
cussed  in  Section  2  below,  but  the  liquid  is  taken  to  have  a  zero-shear- 
rate  viscosity  fiQ.  The  two  axes  of  rotation  intersect  at  the  point  0  on 
the  center  line  of  the  cylinder. 

This  report  is  organized  in  the  following  way.  In  Section  2  we  discuss 
properties  of  non-Newtonian  fluids  that  could  be  relevant  to  the  present 
study,  and  describe  a  number  of  possible  theoretical  models  that  can  be 
used.  The  general  equations  that  govern  the  motion  of  the  liquid,  and 
the  associated  boundary  conditions,  are  developed  in  Section  3,  where 
appropriate  nondimensionalizations  ore  also  Introduced.  Then  in  Section 
U  we  derive  the  expressions  for  the  liquid- induced  moments  acting  on  the 
cylinder,  and  infer  a  general  relationship  between  two  of  the  components 
of  the  moment. 

The  work  performed  in  this  project  has  proceeded  along  two  distinct  but 
related  paths.  On  the  one  hand  we  have  developed  an  analytical  treatment 
of  the  problem  in  the  case  of  a  cylinder  of  infinite  length,  following 
the  procedure  of  Herbert  (1985)  for  the  Newtonian  liquid  case.  The  infi¬ 
nite  length  model  has  validity  when  the  aspect  ratio  of  the  cylinder  is 
large,  which  is  frequently  the  case  in  flight  and  laboratory  tests,  and 
can  predict  the  flow  field  except  close  to  the  cylinder  ends.  A  reason¬ 
able  estimate  of  the  d^spin  moment  can  also  be  obtain*^  ad-ie  procedure 
and  results  are  outlined  in  Sections  5  and  6  respectively.  The  method  of 
analysis  utilizes  the  fact  that  the  nutation  angle  9  is  generally  quite 
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small  (|tf|  S  20°)  ,  and  that  the  ratio  of  angular  velocities  (Q/w)  is 
also  quite  small,  These  facts  enable  the  solution  of  the  problem  to  be 
found  by  analytic  perturbation  techniques,  as  was  done  by  Herbert  (1985) 
for  the  Newtonian  liquid  case,  In  performing  this  analysis  we  were  able 
to  assess  the  contribution  to  the  despin  moment  arising  from  the 
liquid's  viscoelasticity,  and  to  show  how  this  contribution,  from  a 
theoretical  point  of  view,  depends  on  the  constitutive  model  used  to 
represent  the  non-Newtonian  behavior. 

In  parallel  with  the  perturbation  analysis  we  have  performed  a  numerical 
simulation  of  the  flow  in  the  cylinder  shown  in  Figure  1,  and  have  com¬ 
puted  the  moments.  The  computations  were  done  using  a  task-oriented  ver¬ 
sion  of  the  finite-element  code  FIDAP;  this  code  and  the  methodology  of 
the  simulations  are  described  in  Section  7.  The  properties  of  the  non- 
Newtonian  fluids  to  be  simulated  were  deduced  from  data  supplied  by 
CRDC,  as  were  the  parameter  ranges  of  interest.  The  results  of  the 
computations  are  presented  in  Section  8. 

The  results  obtained  are  reviewed  and  discussed  in  Section  9. 
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2.  NON -NEWTONIAN  LIQUIDS 

In  attempting  to  model  the  flows  of  non-Newtonian  liquids  the  funda¬ 
mental  difficulty  always  is  the  form  of  the  relationship  between  the 
stress  tensor  and  the  deformation-rate  (strain-rate)  tensor,  the  so- 
called  constitutive  relation.  No  constitutive  relation  has  been  found, 
either  empirically  or  from  theory,  that  can  claim  to  be  appropriate  for 
all  flows  of  all  non-Newtonian  liquids.  In  fact,  it  is  generally  agreed 
that  there  is  no  such  general  constitutive  relation  or,  if  there  is,  it 
is  so  complex  that  it  cannot  be  used  for  the  solution  of  problems  and 
cannot  be  determined  empirically. 

For  modeling  purposes,  therefore,  it  follows  that  one  is  confronted  with 
the  more  restricted  task  of  determining  a  constitutive  relation  that 
applies  to  a  limited  class  of  flows  and/or  a  limited  class  of  fluids  If 
one  asks  whether  it  is  possible  to  find  a  constitutive  relation  that 
describes  all  possible  flows  of  one  given  non-Newtonian  liquid,  the 
answer,  based  on  experience,  is  again  negative.  It  turns  out,  on  the 
other  hand,  that  there  are  much  better  prospects  of  finding  a  constitu¬ 
tive  relation  to  describe  a  limited  class  of  flows  for  a  reasonably  wide 
range  of  fluids. 

The  consequence  of  this  fact  is  that  in  studying  a  particular  flow  it  is 
of  primary  importance  to  choose  a  constitutive  relation  that  is  appro¬ 
priate  to  that  flow.  One  aspect  of  this  approach,  however,  needs  to  be 
carefully  monitored.  As  will  be  seen  below,  constitutive  relations 
always  contain  parameters  that  measure  the  viscoelastic  properties  of  a 
liquid.  These  parameters  have  to  be  determined  empirically,  but  the 
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de termination  has  to  be  from  experiments  on  the  same  class  of  flows  for 
which  the  model  is  being  constructed,  since  otherwise  they  may  not  be 
relevant  to  the  particular  constitutive  relation  being  applied. 


For  example,  it  is  well  known  that  many  polymer  solutions  exhibit  shear 
thinning  in  plane,  unidirectional  shear  flows;  for  such  flows,  and  when 
no  other  fluid  behavior  is  of  interest,  it  is  often  sufficient  to  model 
shear  thinning  by  introducing  a  shear-rate  dependent  viscosity  function 


M  -  /i(7> 


(2.1) 


where  7  is  the  magnitude  of  the  rate  of  shear.  From  equation  (2.1)  one 


deduces  a  constitutive  law 


£  -  ^(7)2 


(2.2) 


which  is  an  explicit  relation  between  the  stress  tensor  r_  and  the  rate 
of  strain  tensor  j.  Various  functional  forms  of  equation  (2.1)  have  been 
proposed,  mostly  involving  algebraic  dependence  on  7;  a  common  form  is 


the  Carreau  model 


„  f 

**o  -  V 


1  +  (>7 )' 


(2.3) 


where  is  zero-shear-rate  viscosity,  is  inf inite - shear-ra  :e  viscos¬ 
ity,  X  is  a  time  constant  and  a  a  power -law  index. 


While  the  constitutive  relation  implied  by  equation  (2.2)  can  of'  ■  n  be 
fitted  quite  well  to  experimental  data  on  shear  thinning  in  sir-'Ie 
shearing  flows,  it  is  incapable  of  predicting  other  common  viscoelastic 
effects,  in  particular  normal  stress  difference  ar.d  fluid  elasticity. 


Wv-  .^jVudVV'  '•-'j*'  .V  JV  t  s  '  /  T  ^  /-  ■: 
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For  these  and  others  it  is  usually  necessary  to  go  to  a  more  sophisti¬ 
cated  model.  Thus,  for  example,  the  constitutive  relation  known  as  th. 
second-order  fluid  model,  typically  written  in  the  form 

L  “  Qi2  +  a:Di  "  Qiii'2  •  (2.4) 

where  £,  are  the  appropriate  tensors,  c^,  and  c*n  are  constants,  and 
D  represents  a  convected  or  Jaumann  derivative,  predicts  no  shear  thinn¬ 
ing  but  does  predict  a  normal  stress  difference  and  a  complex  viscosity 
in  relaxation.  Again,  it  should  be  emphasized,  models  such  as  (2.4)  can 
bn  made  to  fit  experimental  observations  only  in  a  limited  class  of 
simple  flow  configurations.  An  improved  version  of  equation  (2.4)  is  the 
CEF  model  (Criminale  et  al ,  1958)  in  which  the  constants  a.,  a2,  an  are 
ropiacod  by  (empirically  determined)  functions  of  shear  rate. 

A  class  of  more  widely  acceptable  models  are  the  differential  relations 
of  Che  form 

r  XDr  -  2hq(2  +  eAPjr)  (2.5) 

whir.'/  contain  a  relaxation  time  constant  ,  a  zero-shear  rate  viscosity 
/»0  and  a  retardation  time  constant  e\.  The  derivative  D  includes  a  para¬ 
meter  that  changes  the  character  of  the  model;  the  details  will  be  dis¬ 
cussed  further  below.  This  model  (2.5)  is  known  to  be  applicable  over  a 
wider  class  of  flows  than  any  of  the  others  mem i  •'■•.e  ’'-.eviousiy,  but  it. 
is  also  known  to  be  limited,  and  even  unsatisfi  ,  "tin  applied  to 

very  complicated  flows. 
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In  the  present  project  we  are  dealing  with  a  fully  three-dimensional 
flow  of  a  viscoelastic  liquid.  For  three-dimensional  flows  there  is  vir¬ 
tually  no  information  available,  either  experimental  or  theoretical,  to 
determine  whether  one  or  any  of  the  known  constitutive  relations  is  ade¬ 
quate  to  describe  and  predict  the  flow  dynamics.  This  is  because  it  is 
not  known  whether  shear  thinning,  normal  stress  differences,  elongation- 
al  viscosity,  stress  relaxation,  or  some  combination  of  all  or  some  of 
these,  plays  the  dominant  role  in  any  given  three-dimensional  flow.  On 
this  account  a  major  thrust  of  this  project  has  been  the  testing  and 
categorizing  of  constitutive  relations,  with  a  view  to  eventual  calibra¬ 
tion  against  experiments,  to  ascertain  what  are  the  most  important  char¬ 
acteristics  affecting  the  flow  i  i  a  spinning  and  nutating  cylinder.  The 
object  of  this  analysis  was  to  arrive  at  a  sufficiently  accurate  predic¬ 
tive  constitutive  relation  that  can  be  used  for  large-scale  simulation. 
In  Sections  5  and  6  of  this  work  we  describe  the  results  of  testing 
models  such  as  those  mentioned  explicitly  above,  and  draw  some  relevant 
conclusions . 

The  two  principal  (but  by  no  means  only)  measures  of  the  degree  of  "non- 
Newtonianness"  of  a  liquid  are  the  zero  -  shear- rate  viscosity  /uQ  and  a 
time  constant,  denoted  A.  The  addition  of  certain  polymers  to  a  Newto¬ 
nian  liquid  has  the  effect  of  increasing  the  viscosity,  sometimes  by 
orders  of  magnitude,  so  that  refers  to  the  viscosity  of  the  new,  non- 
Newtonian  liquid  and  is  by  no  means  the  same  as  that  of  the  solvent.  For 
the  case  of  liquid- filled  shells  at  high  (Newtonian)  viscosity,  this 
means  that  the  new  polymer  solution  must  be  regarded  as  having  an  even 


higher  viscosity. 
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The  time  constant  A  is  a  measure  of  the  degree  of  shear  thinning  under 
deformation,  as  well  as  of  normal  stress  and  stress  relaxation  depending 
on  the  type  of  flow  under  consideration.  To  estimate  how  important  any 
or  all  of  these  effects  are  in  any  given  situation,  one  should  compare 
the  magnitude  of  A  (which  has  the  dimensions  of  time)  with  some  other 
typical  time  constant  of  the  flow.  In  the  present  work  an  obvious  candi¬ 
date  for  the  latter  is  the  period  of  the  spin,  in  which  case  we  can  de¬ 
fine  a  dimensionless  parameter  We,  called  the  Weissenberg  number,  by  the 
formula 

We  -  Aw  .  (2.6) 

Elastic  effects  become  important  when  the  elastic  time  constant  is  large 
compared  with  the  dynamic  time  constant;  this  may  well  be  the  case  in 
the  present  situation  since  the  angular  velocity  of  spin  is  quite  large 
(4000-6000  r.p.ra.). 

For  estimating  the  dynamical  behavior  of  the  liquid  the  appropriate 
dimensionless  parameter  is  the  Reynolds  number 


Combining  equations  (2.6)  and  (2.7)  we  can  define  an  alternative  elastic 
paramet'  r,  called  the  Deborah  number,  by  the  formula 


De  -  We/Re  - 


(2.8) 
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which  has  the  virtue  that  it  is  independent  of  spin  rate,  and  represents 
the  ratio  of  the  elastic  time  constant  to  a  typical  diffusion  time 
constant , 

The  parameters  Re,  We,  De  will  play  an  important  role  in  the  subsequent 
discussion. 
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3 .  FORMULATION 


As  pointed  out  by  Vaughn  et  al  (1983)  and  Herbert  (1985) ,  it  ts  conven¬ 
ient  to  write  the  governing  equations  in  a  nutating  frame  of  reference, 
the  so-called  aeroballistic  frame.  In  this  coordinate  system,  the  z-axis 
coincides  with  the  axis  of  the  cylinder,  the  x-axis  lies  in  the  plane 
containing  the  angular  velocity  vectors  w  and  Q,  and  the  y-axis  is  per¬ 
pendicular  to  this  plane.  This  constitutes  a  right-handed  cartesian  sys¬ 
tem,  with  the  origin  taken  at  the  center  of  the  cylinder. 


It  is  also  convenient  to  use  dimensionless  forms  of  the  governing  equa¬ 
tions.  Nondimensionalization  is  achieved  by  scaling  lengths  with  respect 
to  the  radius  a,  time  with  respect  to  the  spin  time  constant  w'1,  linear 
velocities  with  respect  to  the  speed  aw,  the  angular  velocity  vector  Q 
with  respect  to  its  magnitude  0,  the  stress  and  deformation  rate  tensors 
with  respect  to  and  w  respectively,  and  the  pressure  with  respect  to 
pa2ul.  It  is  then  easy  to  show  that  the  governing  equations  for  momentum 
and  mass  conservation  with  respect  to  the  aeroballistic  frame  are 


7  1 

yfc  +  v.  Vy  +  2 rjQ  X  v  +  rjiQ  X(Q  x  r)  -  -Vp  +  —  V-  r 


V-v  -  0  . 


(3.1) 


(3.2) 


The  parameters  appearing  here  are  the  Reynolds  number,  as  defined  by 


formula  (2.7),  and  the  spin  ratio  rj 


t)  -  n/w  . 


(3.3) 


The  dimensionless  angular  velocity  vector  Q  with  respect  to  the  carte¬ 
sian  aeroballistic  frame  is 

Q  -  -ol  +  J~az  k  (3.4) 


where 


a  -  sin  9 


(3.5) 


is  the  third  crucial  parameter  of  the  problem.  The  boundary  conditions 
are  that  the  normal  velocity  component  is  zero  at  every  rigid  boundary 
and  that  the  tangential  velocity  at  a  boundary  is  equal  to  the  rigid- 
body  velocity  of  that  boundary.  On  physical  grounds  we  require  also  that 
the  velocity  is  bounded  on  the  cylinder  axis. 

To  complete  the  specification  of  the  problem  it  is  necessary  to  stipu¬ 
late  the  constitutive  relation  between  stress  and  strain  rate.  We  shall 
consider  a  number  of  typical,  representative  relations  all  of  which,  it 
can  be  shown,  are  invariant  with  respect  to  the  transformation  from  the 
inertial  frame  to  the  aeroballistic  frame. 


The  dimensionless  form  of  the  constitutive  relation  corresponding  to  the 
power-law  Carreau  model  (2.3)  is 


2i 


[«  +  (1-0(1  +  We  r 


(3.6) 


where  We  is  defined  by  equation  (2.6),  t  -  /b0  and  is  often  taken  to 
be  zero,  and  a  is  the  power-law  index.  The  strain-rate  tensor  in  (3.6) 
is  defined  by 


2  “  ^(Vv  +  Vv1) 


(3.7) 


and  7  is  its  magnitude. 


^•v*  v « *s ts  ■a  'ATVTi  ww  *  v*  v*y  •ir‘^x*x*xmy^jrrjny*¥*y*jrix”xvx*y?<r,jr>y'*'my?}'>  ??  v 
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The  dimensionless  form  of  the  second-order  fluid  relation  (2.4)  is 

r  -  2(j  -  We  +  2 *  (3.8) 

where  We  is  the  Weissenberg  number  again  and  k  is  a  constant.  The  deriv¬ 
ative  D  appearing  in  (3.8)  is  defined  for  a  second-order  tensor  p  by  the 
formula , 

Dp  -  pt  +  v-Vp  +  wp  -  p’<*>  -  a(- y-p  +  p--y)  ,  (3.9) 

where  w  is  the  vorticity  tensor 

tt  -  |(?Y  -  VvT)  (3.10) 

and  a  is  a  constant  that  can  take  values  between  -1  and  +1.  The  choices 
a  -  -1,0, +1  reduce  equation  (3.9)  to  the  lower  convected,  corotational 
and  upper  convected  derivative  respectively.  It  should  be  noted  that 
(3.8)  predicts  no  shear  thinning  in  simple  shear  flow,  but  does  predict 
constant  first  and  second  normal  stress  differences.  The  Weissenberg 
number  that  appears  in  (3.8),  consequently,  is  a  measure  of  the  first 
normal  stress  difference  and  therefore  may  need  to  be  interpreted  dif¬ 
ferently  from  the  corresponding  parameter  in  (3.6),  where  it  is  a  mea¬ 
sure  of  shear  thinning.  The  parameter  k  in  (3.8)  is  associated  with  the 
second  normal  stress  difference. 

A  generalization  of  (3.8)  is  the  Criminale - Ericksrn- Filbey  model,  namely 
r  -  2(4>1(-r)2  ~  We4>2(-y  )Djr  +  2x^(7)^  ■  -  (3.11) 

where  <t  $>2,  <t>3  are  generally  •  npiricaily  determined  functions  of  y. 


-.’  VI.  ■.-  «ii  V  V-  -.1  -.-  V  .,-  v- 


V\  *,  -\  -.A  VA  WH.’,  i.\  ■t  -V  »fc^« 


In  particular,  a  suitably  chosen  form  of  can  describe  shear  thinning 
effects ,  while  the  dependence  of  $2  and  $>3  on  7  can  describe  the  varia¬ 
tion  of  first  and  second  normal  stress  differences  with  shear  rate. 

A  rather  general  differential  model  which  we  propose  to  study  is 


r  +  We  Dr  —  2(jr  +  <We  Djr) 


(3.12) 


where  the  operator  D  is  defined  by  (3.9).  This  model  predicts  shear 
thinning  when  -1  <  a  <  1  and  a  shear -rate  dependent  first  normal  stress 
difference.  It  is  a  modified  version  of  the  Oldroyd  (1958)  8-constant 
model .  The  Weissenberg  number  in  this  model  is  a  measure  of  both  shear 
thinning  and  normal  stress,  while  t  -  /i  //i0,  as  in  (3.6). 


It  is  convenient  to  introduce  cylindrical  polar  coordinates  (r,^,z), 


defined  by 


x  -  rcos^,  y  -  rsir<£,  z  -  z 


(3.13) 


In  that  case  the  boundary  conditions  can  be  written  simply  as 


v  —  4>  on  r  -  1 


v  -  r<f>  on  z  —  ±  b 


(3.14) 


v  bounded  as  r  -*  0 


where  b  -  c/a  is  the  aspect  ratio,  while  the  angular  velocity  vector 
(3.4)  takes  the  form 


/  ? 

0  -  -ocos<(>  r  +  osin<£  <f>  +  /l -a  k 


(3.15) 


This  completes  the  formulation  of  the  problem. 


T ^  «-  V’  **-  ,w%f V »  V* 
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4 .  MOMENTS 

In  this  Section  we  discuss  some  issues  concerning  the  moment  acting  on 
the  cylinder  shown  in  Figure  1,  and  derive  an  important  relationship 
between  two  components  of  the  moment  vector.  In  the  latter  part  of  the 
Section  we  show  what  modifications  are  required  when  the  problem  of  the 
infinitely  long  cylinder  is  being  considered. 

The  total  moment  acting  on  the  cylinder  is  denoted  M^.  Typical  units 
are  gmcm2- sec’2,  and  in  conformity  with  the  scalings  introduced  in  Sec¬ 
tion  3  we  define  the  dimensionless  total  moment  Mj  by  the  formula 

fcf*  -  paV  Mt  (4.1) 

where  p  is  fluid  density,  a  is  cylinder  radius  and  w  is  the  spin  rate. 

An  expression  for  the  moment  can  be  written  down  in  terms  of  pressure 
and  stress  fields,  namely 

Mt  -  Jr  x(-Vp  +  ^  V-r)dV  (4.2) 

or,  equivalently, 

Mt  -  J(-pr  x  dS  +  £  x  r  ■  dS)  (4.3) 

where  V  denoted  the  volume  of  liquid  and  5  the  bounding  surface.  With 
the  aid  of  (3.1)  it  is  also  possible  to  w’-ite  the  moment  in  the  form 

M  =  j-  J(r  x  v).iV  +  J(r  x  v)  (v  dS)  +  2r?Jr  x(Q  x  v )dV 

V  S  V 


+  r7ZJn  x(Q  X  r)dV 


(4.4) 
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with  respect  to  the  aeroballistic  coordinate  system.  This  formula 
expresses  the  fact  that  the  moment  is  equal  to  the  rate  of  change  of 
angular  momentum, 

The  fourth  terra  on  the  right-hand  side  of  (4.4)  represents  the  moment 
due  to  solid-body  motion  and  is  Independent  of  the  velocity  field.  It  is 
therefore  convenient  to  remove  this  term  and  to  introduce  the  quantity 
M,  representing  the  (dimensionless)  moment  Induced  by  the  liquid  pay- 
load,  defined  by 

M  -  —  J(r  x  v)dV  +  J(r  x  y)(ydS)  +  2qjr  x(Q  x  y)dV  (4.5) 

dt  v  s  v 

The  dimensional  equivalent  of  this  is  given  by 

il*  -  pasu2  £1  .  (4.6) 

The  solution  being  sought  is  time- independent  with  respect  to  the 
aeroballistic  reference  frame;  this  implies  that  the  first  term  on  the 
right-hand  side  of  (4.5)  is  identically  zero.  Moreover,  the  boundary 
conditions  (3.14)  imply  that  v-dS  -  0  on  the  entire  bounding  surface.  It 
follows,  therefore,  that  the  effective  form  of  (4.5)  is 

M  -  2r;Jr  x(Q  x  v )dV  (4.7) 

v 

which  is,  of  course,  the  contribution  from  the  coriolis  force. 


We  refer  the  moment  to  a  cartesian  coordinate  system  in  the  aeroball¬ 
istic  frame;  this  is  the  system  described  at  the  beginning  of  Section  3. 
The  components  of  r  are  (x,y,z),  while  the  components  of  M  are  denoted 


K, W  v~»  V*' •'ft-'. 


Mx,  rty,  Mt.  In  this  terminology  Mz  represents  the  despin  (roll)  moment, 
the  yaw  (side)  moment  and  My  the  pitch  moment.  For  the  velocity  vec¬ 
tor  v  it  is  convenient  to  write 


y  -  (-y  +  u)i  +  (x  +  v)j  +  wk  , 


(A. 8) 


so  that  (u,v,w)  are  the  components  of  the  deviation  from  rigid-body 
rotation.  It  should  be  noted  that  the  boundary  conditions  (3.14)  imply 
the  vanishing  of  u,  v  and  w  on  all  solid  boundaries.  The  components  of 
the  angular  velocity  vector  Q  are  given  by  (3.4). 

Substituting  all  the  various  components  into  (4.7)  we  obtain 
Mx  -  -2r;J*(yvsin?  +  zvsin?  +  zucosfi)dV 


(4.9) 


My  -  rjanb  +  2fjJ(xvsin?  -  zvcos?)dV 


Mz  -  2»jJ"(xwsin?  +  xucosff  +  yvcos8)dV 


where  8  is  the  coning  angle  and  b  -  c/a  is  the  aspect  ratio.  It  is  note¬ 
worthy  that  the  rigid-body  motion  of  the  liquid  contributes  only  to  the 
pitch  moment  N . 

We  now  establish  an  exact  relationship  between  the  components  and  M  . 
To  do  this  it  is  convenient  to  write  the  first  and  third  of  equations 
(4.9)  in  the  forms 


Mx  -  -2i7(sin0  I1  +  sin?  I2  +  cos?  I3) 


-  2rj(sin?  Kx  +  cos?  K2  +  cos 9  K3) 


(4.10) 


> 1  J - to* -  tfcj 1  '-r>  to>  *'»  W**  V ?K 


where  the  definitions  of  the  I's  and  K's  are  self-evident.  Consider 
first 

l1  -  Jyv  dV  -  JXfyvdxdydz  -  \  JJ‘Jvd(y2)dxdz  . 

v 

Integrating  by  parts  and  using  the  fact  that  v  -  0  on  the  boundaries,  we 
have 


by  continuity.  Now  ve  can  write 

lr  -  I JJV2  S  ^  dxdydz  +  I JJV2  J  fi  dzdydx 

and  the  interior  integrals  both  vanish  since  u  -  0  and  v  -  0  on  the 
boundaries.  Hence  Ix  -  0. 

By  exactly  similar  procedures  we  can  show  that  I2,  and  K3  also  van¬ 
ish.  Thus  we  have  the  result 

lx  -  I2  "  K2  -  K3  -  0  .  (4.11) 

This  implies  that  equations  (4.10)  reduce  to 

Hx - 2rjcosfl  I3  ,  Mt  -  2f?sin0  Kj  .  (4.12) 

Next  consider 

i3  -  <*V  -  JJz  Ju dxdydz  -  -Jzx~  dV 

v  u 


on  integrating  by  parts  and  using  the  condition  u  -  0  on  the  boundaries. 


-19- 


The  mass  of  liquid  inside  the  control  volume  is 

-  2w azcp  (4.15) 

so  that  we  can  conveniently  define  a  liquid-payload- induced  moment  per 
unit  mass  by  the  formula 

m*  -  M*/^  ;  (4.16) 

this  has  units  cm2- sec’2.  A  dimensionless  moment  per  unit  mass,  denoted 
m,  can  be  obtained  from  this  by  writing 

m*  —  azu2  m  (4.17) 

and  it  is  then  easy  to  show  that  m  and  U  are  related  by  the  formula 

m  -  H/(2 nb)  .  (4.18) 

The  quantities  ffi*  and  m  are  meaningful  for  both  a  finite  cylinder  and  an 
infinite  cylinder,  and  therefore  can  be  used  when  results  for  the  two 
problems  are  being  compared. 

The  second  modification  is  necessitated  by  the  fact  that  the  flat  ends 
of  the  control  volume  are  not  rigid  boundaries,  and  therefore  it  is  no 
longer  the  case  that  v ■  dS  -  0  on  all  boundaries.  It  follows  that  the 
effective  form  of  (4.5)  is  no  longer  (4.7),  but  rather 

M  =  J(r  x  v)  (v •  dS)  +  2r?Jr  x(Q  x  v)dl'  ,  (4.19) 

s  v 

where  the  first  term  on  the  right-hand  side  of  this  equation  may  have 
contributions  from  the  flat  ends.  Precisely  what  these  contributions  are 


will  become  appatcnC  in  the  next  Section  where  the  calculations  for  the 

1  i 

infinitely  long  cylinder  are  pre -anted. 

Finally,  we  do  not  expect  the  relationship  (4.14)  to  hold  in  the  infin¬ 
ite  cylinder  case.  Not  or./y  are  the  moments  affected  by  the  additional 
term  in  (4,19\  but  also  the  identities  (4.11)  and  (4.13)  cannot  be 
established.  Thit  ir.  because  the  proofs  of  these  depend  on  the  fact  that 
the  components  are  all  zero  on  all  boundaries;  but  this  is  no 

longer  true  on  the  flat  ends  of  the  control  volume.  Consequently  the 
proofs  break  down  and  the  results  do  not  apply.  In  the  next  Section 
shall  discuss  the  relationship  between  and  as  it  emerges  from  our 


solution. 


jnurMr»jr*jf*s  *w*  vr^mx  *jr Kmu j»  >i  ^y  >*  ax  r  v  -v»  ax  -m  r.i  r.*  ^  *v*  a*  ,m»  a*  -v*  a*  .-a* . 
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5.  INFINITELY  LONG  CYLINDER 


In  this  Section  and  the  next  we  shall  be  considering  the  special,  and 
obviously  unrealistic,  case  of  an  infinitely  long  cylinder.  This  problem 
was  analyzed  in  considerable  detail  by  Herbert  (1985)  for  the  Newtonian 
liquid;  we  shall  adapt  Herbert's  procedures  to  the  non-Newtonian  liquid, 
and  shall  indicate  similarities  and  dissimilarities  as  appropriate. 


There  are  two  reasons  for  pursuing  the  inf inite -cylinder  problem.  The 
first  is  that  it  is  possible  to  obtain  an  approximate  analytic  solution 
in  the  case  that  the  spin  ratio  (O/w)  is  small,  which  does  indeed  apply 
in  most  situations  of  practical  interest.  The  second  is  that  the  results 
obtained  should  be  relevant  to  the  case  of  a  cylinder  which  is  finite 
but  of  large  aspect  ratio.  As  indicated  in  the  previous  Section, 
however,  and  shown  further  below,  some  reservations  need  to  be  made  in 
interpreting  inf  inite -cylinder  results  for  the  moments  when  one  has  the 
finite  cylinder  in  mind. 

The  governing  equations  are  (3.1) -(3. 2),  together  with  one  of  the  con¬ 
stitutive  relations  listed  ir  ’action  3.  The  boundary  conditions  (3.14) 
are  modified,  however,  for  infinite  aspect  ratio  to  read 


A 

v  -  <f>  on  r  -  1 

v  bounded  as  r  -»  0 


(5.1) 


It  is  assumed  that  the  spin  ratio  q  is  small,  and  a  solution  is  devel¬ 
oped  through  a  regular  perturbation  procedure. 


V  *j  A:  /.  V.  WV  -A.  , 


We  expand  the  field  quantities  appearing  in  the  governing  equations  in 
power  series  of  the  following  forms, 


v  -  y<0)  +  +  r)iwlZ)  + 


r  -  r<0)  +  +  r?V2>  + 


p  -  p<0)  +  f,pa)  +  r,2p<2>  + 


(5.2) 


Associated  with  these  will  be  similar  expansions  of  the  strain-rate  and 
vorticity  tensors, 


2  -  2<0>  +  V21"  +  r,2Zi2)  +  ••• 


(5.3) 


,,  .  ,,<°>  4.  ,C1 )  ,  _2  (2)  . 

CJ  “  (J  T  TJ(j>  ~r  T)  (j)  + 


These  representations  are  substituted  into  the  governing  equations  (3.1) 
-(3.2),  the  boundary  conditions  (5.1),  and  the  appropriate  constitutive 
equations,  and  the  terms  corresponding  to  like  powers  of  rj  are  equated. 
This  results  in  a  sequence  of  boundary-value  problems  which  are  in 
principle  amenable  to  analytic  solution. 


It  is  easy  to  show  that  the  leading-order  solution,  corresponding  to 
*1  -  0,  is  simply  the  rigid-body  motion  of  the  liquid,  and  this  is  inde¬ 
pendent  of  the  particular  constitutive  equation  chosen.  In  fact  we  have 
that 

v(0)  -  4,  r(0)  -  0,  p(0>  .|r!, 

~  (5.4) 


,<0) 


r  <f>  -  d  r 


It  should  be  noted  that  the  strain-rate  tensor  is  zero,  while  the 
vorticity  tensor  has  components  appropriate  to  rigid-body  rotation. 


At  the  next  order,  corresponding  to  0(r;),  equations  (3.1) -(3. 2)  become 


vf”  +  v"’  +  (v(1)  -V)r£  +  2n  x  rj  -  -7p(1)  +  V.r(1> 

— t  —  *  —  —  —  r  Re  - 

V.ya>  -  0 

while  the  boundary  conditions  (4.1)  become 

v(1)  -  0  on  r  ■  1 
va)  bounded  as  r  -*  0  . 


(5.5) 

(5.6) 

(5.7) 


The  appropriate  approximation  for  each  of  the  constitutive  relations 
discussed  in  the  preceding  Section  can  easily  be  written  down.  From  the 
Carreau  model  (3.6)  we  obtain 


(5.8) 


which  is  the  same  a3  for  a  Newtonian  liquid.  This  is  due  to  the  fact 
that  7  appears  in  Che  denominator  of  (3.6)  in  quadratic  form,  and  there¬ 
fore  the  nonlinearity  makes  no  contribution  at  this  order. 

For  the  second-order  fluid  model  (3,8)  we  obtain 


r(1>  _  22(1)  -  2We(2‘n  +  i}1’  +  w(t»-  2U’ 


2ll>-  “‘0’' 


(5.9) 


which,  it  should  be  noted,  does  not  contain  either  of  the  parameters  a 
or  k.  The  form  of  the  CEF  equation  (3.11)  depends  on  the  functions  C>  , 
$2 ,  <X>3  but  if,  as  is  usually  the  case,  these  are  quadratic  in  7,  then 
the  reduced  form  of  the  CEF  equation  is  also  (3.9),  the  same  as  that  of 
the  second- order  fluid  equation. 
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The  differential  model  (3.12)  becomes 


r<l>  +  We^1’  +  r^15  +  a>(0)-r<1)  -  r(1).w<0>} 


(5.10) 


-  2ia)  +  2€We(2‘1)  +  +  «"”•  j14’  -  2llJ.  w'1”) 


(1)  ,  •  (1)  ,  ,.<0)  H)  fl)  (0), 


It  should  be  observed  that  the  parameter  a,  which  determines  the  extent 
of  shear  thinning  in  the  model,  is  absent  from  (5.10);  in  this  problem 
upper  convicted,  lower  convected  and  corotational  models  all  reduce  to 
the  same  form  in  this  approximation. 


The  problem  determined  by  equations  (5. 5) -(5. 6)  and  one  of  (5.8),  (5.9) 
or  (5.10)  has  a  time  - independent  solution  with  the  velocity  vector  being 
of  the  form 


Y<1>  -  t  0,  0,  v<r,*>  ]  . 


(5.11) 


The  flow  is  purely  axial  and  depends  only  on  the  plane  coordinates  r,4>. 
With  this  flow  field  the  strain-rate  and  stress  tensors  are  respect¬ 
ively, 
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(5.12) 


Equations  (5.5)  and  (5.6)  now  simplify  to 


pu>  -  J\-oz ■  r2 


(5. ' D; 


and  the  equation, 


-  2or  cos*  -  ^  (i(rrl3)r  +  *7^ 


(5.14) 


<&■  ttetmW  »'i  ^»  f«  ^  r-d  *>,  4*  ft4.4V4 
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with  che  boundary  conditions 


w  -  0  on  r  -  1,  w  bounded  at  r  -  0  . 


(5.15) 


In  the  case  of  the  (Newtonian)  constitutive  relation  (5.8)  we  find  that 


f  a  tj 

13  r’ 


T  -  If. 

23  t  <t> 


(5.16) 


so  that  (5.14)  becomes 


v.  - 


2 or  cos^  -  7*w 


(5.17) 


where  72  is  the  two-dimensional  Laplacian  operator  in  the  r,^-plane. 
This  is  the  equation  solved  by  Herbert  (1985). 

In  the  case  of  the  second-order  fluid  relation  (5.9)  we  find  that 


r  13  “  ("  -  We  W*>r  *  r23  “  J<w  -  We  Wt\ 


(5.18) 


which  reduces  (5.14)  to  the  equation 


1  2 

w.  -  2ar  cos^  -  —  V  (w  —  We  w.) 
*  Ke  * 


(5.19) 


In  the  case  of  the  differential  model  (5.10)  we  obtain 


ri3  +  We  A  *  (V  +  «We 


13,  £ 


r23  +  We  f23,*  "  7(W  +  <Ue  V* 


(5.20) 


Although  these  equations  can  be  solved  for  Tn  and  r23,  and  the  solution 
substituted  into  (5.14),  we  shall  refrain  from  doing  so  at  this  point 
since  the  subsequent  development  provides  an  easier  way  of  achieving  the 
desired  result. 
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Following  Herbert  (1985)  we  solve  the  problem  by  writing 


w  -  2a  [f(r)cos<£  +  g(r)sin^) 


-  2a  ■  Real  (h(r)e'1*] 

where  h  -  f  +  ig.  Equation  (5.17)  becomes 

r2h"  +  rh‘  -  <1  -  iRer2)h  -  -Rer3 


Equation  (5.19)  becomes 


r2h"  +  rh' 


i  Rer2 
1+iWe 


)” 


-Rer3 

1+iWe 


(5.21) 


(5.22) 


(5.23) 


In  the  third,  differential  case  we  find  the  solutions  of  (5.20)  to  be 


r 


13 


r 


23 


2a -Real| 
2a-Real[ 


(l-lfWe)h'  c1* 
1  -  IWe 


-i(l-i  fWe)he-1*  1 
r(l  -  iWe)  j 


(5.24) 


whereupon  the  appropriate  equation  is  found  to  be 


r2h"  +  rJl.  _  [i  _  iR'(irW±).£*  lh  „ 

(  1  -  ie We  J  1  -  r<We 


(5.25) 


The  boundary  conditions  are 


h  -  0  at  r  -  1 ,  h  bounded  at  r  -  0 


(5.26) 


We  see  that  equations  (5.22),  (5.23)  and  (5.25)  all  have  the  same  struc¬ 
ture.  In  fact  they  can  all  be  written  in  the  single  unified  form 

z2h"  +  rh'  -  (1  -  lSr2)h  - Sr3  (5.27) 


Q  s'Vs" 


*»•  v  v.**  V  yr  V  «- 


v'< 
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where 


S 


S 


Re 

Re 

1+iWe 


(Newtonian  and  power-law  models) 


(second-order  fluid  model) 


(5.28) 

(5.29) 


s  _  (l-|We)Re  (differential  model)  (5.30) 

i - l ewe 

It  is  noteworthy  that  for  the  second-order- fluid  and  differential  models 
the  parameter  S  has  the  form  of  a  complex  Reynolds  number.  Thus  in  these 
two  cases  it  seems  that  the  non-Newtonian  effects  are  in  a  sense  repre¬ 
sentable  by  a  complex  viscosity;  this  occurs  in  other  flow  configura¬ 
tions  of  viscoelastic  fluids,  for  example  in  time -dependent  flows  where 
stress  relaxation  is  important. 

The  solution  of  (5.27)  which  satisfies  the  boundary  conditions  (5.26)  is 

f  M^)  1 

h-V--i1T5)l  <53L) 

where 

q2  -  -LS  (5.32) 

and  I1  is  the  modified  Bessel  function  of  order  unity.  The  flow  field 
can  be  completely  determined  by  a  simple  computation  of  the  Bessel  func¬ 
tions  . 


Before  proceeding  to  discuss  the  moment  acting  on  the  cylinder,  we  note 
that  approximate  forms  of  (5.31)  for  small  and  large  |s|  are  easily 
found . 
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When  |s|  «  1  we  apply  standard  formulas  for  series  expansions  of  the 
Bessel  functions  to  obtain 

h  -  |(r  -  r3)  +  ff^(2r  -  3r3  +  r5)  +  0 (S3)  .  (5.33) 

In  the  case  of  a  Newtonian  or  power-law  fluid  we  have  that  S-Re ,  so  that 
| S |  «  1  corresponds  simply  to  a  small  Reynolds  number  approximation. 
The  separation  of  (5.33)  into  its  real  and  imaginary  parts  (h  -  f  +  ig ) 
gives 

f  -  |£(r  -  r3)  ,  g  -  ff^(2r  -  3r3  +  r5)  (5.34) 

to  a  leading  approximation,  in  agreement  with  Herbert  (1985). 


For  the  second-order  fluid  (5.29)  the  condition  |s|  «  1  is  achieved 
when  Re  «  1.  It  can  also  be  achieved  for  very  large  We,  but  this  is  an 
anomalous  case  since  the  second-order  fluid  model  is  known  to  be  invalid 
at  large  Weissenberg  number.  The  real  and  imaginary  parts  of  (5.33)  are 
found  to  be 


f  - 


Re 


8(1+We2) 


,  3  x  .  Re2We 

-(r  -  r  )  +  - 


zs  2 


96(1+We  ) 


(2r  -  3r3  +  r5) 


(5.35) 


8  - 


-ReWe 
8 ( 1+Wez) 


(r  -  rJ)  + 


Rez(l-We2) 

192(l+We2)2 


(2r  -  3r3  +  r5) 


this  is  significantly  different  from  (5.34)  in  that  the  component  g  now 
has  a  term  proportional  to  Re,  whereas  there  is  no  such  term  in  equation 


(5.34). 


-29- 


For  the  differential  model  (5.30)  the  condition  |s|  «  1  is  again 
achieved  when  Re  «  1.  Provided  t  h  0  the  condition  Re  «  1  implies  that 
|s|  «  1  for  the  whole  range  of  values  of  We.  When  e  -  0,  however,  it  is 
clearly  necessary  to  impose  the  additional  requirement  WeRe  «  1  to  en¬ 
sure  that  | S |  «  1.  In  this  case  (5.33)  gives 


f  - 


Re(l+£We2)  ,  3,  .  Re2We(l— e) (1+eWe2) 

8(l+£2We2)  96(l+£2We2)2 


(5.36) 


-(1-QReWe 

8(l+ezWe2) 


r3) 


R.^l+.w.y-o-oW)  _  3rS  +  rS 

192(l+£2We2)2 


Again  there  is  a  term  proportional  to  the  Reynolds  number  in  the  compo¬ 
nent  g. 

In  all  three  cases  (5 . 28)  -  ( 5 . 30)  the  limit  |s|  -*  «  is  equivalent  to 
Re  -*  co,  irrespective  of  Weissenberg  number.  The  effects  of  elasticity 
disappear  in  this  limit  and  all  the  cases  are  effectively  Newtonian. 
Asymptotic  results  in  the  limit  Re  -♦  «  have  been  provided  by  Herbert 
(1985) . 

We  proceed  now  to  calculate  the  moment  in  the  present  configuration.  We 
take  a  control  volume  V  with  bounding  surface  S  which  coincides  with  a 
cylinder  of  liquid  of  dimensionless  length  2b  and  dimensionless  radius 
1,  and  whose  center  coincides  with  the  origin  of  the  coordinate  system. 
As  shown  in  Section  4,  equation  (4.19),  the  liquid- induced  moment  acting 
on  this  control  volume  has  the  form 

H  “  Us  +  tlv  (5.37) 
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where 

Us  -  /(£  x  v)(y  dS)  ,  M,  -  2»»fl  x(Q  x  v)dV  .  (5.38) 

s  v 

As  indicated  previously,  the  contribution  Mg  arises  entirely  from  the 
fact  that  the  boundaries  at  the  flat  ends  of  the  control  volume  are  not 
rigid. 

We  shall  compute  Mg  and  My  separately  for  the  case  of  small  spin  ratio. 
First,  we  write 

Ms  -  M^0)  +  +  rj2M^2)  +  •••  (5.39) 

Then,  substituting  (5.2)  and  (5.39),  and  equating  coefficients  of  like 
powers  ofrj,  we  obtain 

M^0)  -  J(r  x  y<0))(y<0).dS)  (5.40) 


Mg11  ■  J(r  x  v(0>)  (vll)  ■  dS)  +  J(r  x  vu>)  (v(0>  •  dS)  (5.41) 


Mg21  -  J(r  x  v(0))  (v<2)  ■  dS)  +  J(r  xvu,)(y  a,dS) 


(5.42) 


+  J(r  x  v<2))  (vl0>  ■  dS) 


for  the  first  three  terms  in  the  series  for  Mg. 


Now  v(0)  is  given  by  (5.4),  and  it  is  obvious  that  v(0)-dS  vanishes 
everywhere  on  S.  Hence  (5.40)  gives 


Mi0)  -  0 


(5.43) 
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Next  consider  (5.41).  For  the  reason  just  stated  the  second  integral  on 
the  right-hand  side  of  (5.41)  is  identically  zero,  while  the  first  inte¬ 
gral  reduces  to  ar.  integral  over  the  two  flat  ends,  with  (5.11)  giving 
ya)-dS  -  ±  wtrdrd4>  on  z  -  ±  b  respectively.  We  then  obtain 

Mg1' - 2bJJ"r2w  £  drd<j>  .  (5.44) 

the  integration  being  over  0  <  r  <  1,  0  <  <j>  <  2n .  We  express  the  moment 
in  terms  of  its  cartesian  components  by  using  the  transformations 


r  -  cos^  i  +  sin^  j 


4>  -  -sin^  i  +  co s<j>  j 


(5.45) 


and  the  representation  (5.21)  for  w.  Then  we  obtain  from  (5.44)  non-zero 
x  and  y  components , 


M*1'  -  -4jrabJr2f(r)dr 


rt'1’ - 4irobJr2g(r)dr  , 


(5.46) 


(5.47) 


and  the  z-component  is  found  to  be  identically  zero. 


It  is  easy  to  show,  finally,  that 


MlZ)  -  0 
— s 


(5.48) 


for  the  following  reasons.  Although  v(2)  has  not  been  calculated  explic¬ 
itly,  it  can  be  demonstrated  with  little  difficulty  that  y(Z)  has  no 
component  in  the  z-direction.  Hence  v'Z)-dS  -  0  on  the  surface,  and  so 
the  first  integral  on  the  right-hand  side  of  (5.42)  vanishes.  The  thiro 
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r, 


i 


integral  vanishes  because  y(0>-dS  -  0  on  the  surface.  The  second  inte¬ 
gral  can  be  shown  to  vanish  because  of  symmetry  considerations. 

Summarizing  these  results  we  have  that 

Ms  =*  Tj  f  -4jrabJ*r2f (r)dr  i  -  4*abJ*r2g(r)dr  1  +  0(r;3)  (5.49) 

i  o  “  o  ~  J 

We  now  turn  to  the  contribution  which  arises  from  the  coriolis 

force.  Writing 

Mv  “  ttj0>  +  »<1>  +  n2f^2)  +  •••  (5.50) 

and  substituting  into  (5.38),  we  obtain 


s, 


<• 


K* 

bl 


-  o 

M^1’  -  2ji  x(Q  x  v(0))dV 
v 

U^2)  -  2ji  x(Q  x  v(U)dW 
v 

Using  (3.15)  and  (5.4)  we  easily  show  that 


(5.51) 


(5.52) 


(5.53) 
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Combining  (5.49)  and  (5.56)  we  obtain 


-  r)  ^-4xabjr2f  (r)drj  i  +  tj^xbo  -  4jrb<7j*r2g(r)drj  j 


(5.57) 


+  rj2- 8xbo2Jr2f  (r)dr  k  +  0(r)3) 


This  can  be  written  in  an  alternative,  more  convenient,  form.  By  inte¬ 
grating  equation  (5.27)  and  using  the  boundary  conditions  (5.26)  we  find 
that 


Jr2h(r)dr 

o 


ib' (1)  .  i 
S  4 


(5.58) 


We  take  the  real  and  imaginary  parts  of  this  expression  and  substitute 
into  (5.57),  which  then  becomes 


M  -  r}  ■  -hixab  Reel 


mi 


+  fj--4xb(7  Imag 


mi 


(5.59) 


+  rj2-8»rb<72  Real^^~i^j  k 


+  0(nJ)  . 


If  the  components  of  M  are  denoted  by  Mx,My,Mz,  we  see  from  (5.59)  that 

Mz  -  -2rjaMx  (5.60) 

which  agrees  with  the  result  obtained  by  Herbert  (1985)  for  the 
|  Newtonian  fluid. 

It  is  important  to  note  that  the  despin  moment  Mz,  which  derives  from 
|  the  moment  of  the  coriolis  force  fly,  can  give  a  reasonable  approximation 

•  to  the  corresponding  quantity  for  a  cylinder  of  finite  length  ( .'1 th 

l 

■  large  aspect  ratio).  The  other  two  components  M  , M  given  by  (5.59)  are 

\ 

i 

! 
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entirely  spurious  as  far  as  the  finite-cylinder  situation  is  concerned, 
since  they  both  contain  contributions  from  surface  integrals  over  the 
flat  ends.  Similarly  the  relation  (5.60),  which  differs  from  (4.14),  is 
also  spurious. 

From  (5.59),  therefore,  we  note  for  future  reference  that  the  despin 
moment  is  given  by  the  expression 

Mt  -  r?2' 8?:bg2  Real^-— (5.61) 


to  leading  order. 

In  the  next  Section  we  give  numerical  results  based  on  the  calculations 
presented  above . 


.  •  A  t  .  f  ,.  V*. ,  «V 
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6.  INFINITELY  LONG  CYLINDER:  RESULTS 


The  purpose  of  the  calculations  reported  in  this  Section  was  to  demon¬ 
strate  how  the  inclusion  of  non-Newtonian  effects  modified  the  flow 
field  and  the  despin  moment  as  computed  by  Herbert  (1985)  .  We  were 
interested  primarily  in  sensitivity  to  departure  from  Newtonian  behavior 
rather  than  determining  even  an  approximation  to  the  solution  that  was 
applicable  to  the  case  of  a  closed  cylinder  of  finite  length.  The  most 
that  can  be  expected  of  the  present  calculation  would  be  some  qualita¬ 
tive  insight  into  the  consequences  of  including  viscoelasticity  in  the 
model . 

The  results  herein  are  based  on  equation  (5.31),  which  gives  the  compo¬ 
nents  f  and  g  of  the  axial  velocity  in  complex  form,  and  on  equation 
(5.61),  which  gives  the  despin  moment  at  leading  order.  The  computation 
of  the  Bessel  functions  involved  is  straight-forward,  and  we  have  used  a 
combination  of  the  series  and  asymptotic  formulas  for  these  functions  to 
cover  the  required  parameter  range. 

We  note  first  that,  according  to  equation  (5.28),  the  power-law  model 
predicts  the  same  velocity  field  and  despin  moment  as  does  the  Newtonian 
model.  It  should  be  emphasized  that  this  prediction  is  entirely  a  conse¬ 
quence  of  the  linearizing  expansions  (5.2)  and  (5.3),  and  would  fail  to 
hold  when  these  expansions  break  down.  This  question  will  be  discussed 
in  more  detail  subsequently. 

For  the  differential  model,  with  "complex  viscosity"  S  given  by  formula 
(5.30),  we  show  in  Figures  2  and  3  respectively  the  variation  of  f  and  g 
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with  radius  at  a  Reynolds  number  Re  -  15.  Included  in  these  Figures  are 
the  appropriate  Newtonian  curves  (We  -  0)  which  are  identical  with  those 
presented  by  Herbert  (1985) .  The  three  non-Newtonian  cases  shown  in 
these  Figures  correspond  to  various  values  of  Deborah  number  De  and 
retardation  parameter  e,  namely,  De  -  0.1,  e  -  0.1;  De  -  0.2,  e  -  0.1; 
De  -  0.2,  <  -  0.2,  respectively.  It  is  clear  that  the  velocity  fields 
are  considerably  distorted  by  comparison  with  the  Newtonian  case,  par¬ 
ticularly  for  the  component  f  which,  from  (5.21),  corresponds  to  the 
flow  in  the  plane  <f>  -  0.  There  is  here  a  backflow  near  the  cylinder  cen¬ 
ter  and  the  appearance  of  a  boundary  layer  near  the  wall. 

Figures  4  and  5  show  the  variation  of  the  same  quantities,  f  and  g,  with 
radius  for  the  same  parameter  values,  but  at  Reynolds  number  Re  ««  50. 
Again  the  distortions  In  the  velocity  fields  for  the  viscoelastic  cases 
are  self-evident.  In  a  very  broad  sense  the  non-Newtonian  velocity 
fields  at  Re  -  15  are  comparable  with  the  Newtonian  velocity  fields  at 
Re  -  50,  which  is  consistent  with  the  shear- thinning  nature  of  the 
viscoelastic  model. 

The  despin  moment  as  a  function  of  Reynolds  number  for  this  model  is 
shown  in  Figure  6.  The  quantity  actually  plotted  in  this  diagram  is,  for 
convenience,  denoted  H  and  defined  by 


i 


M 


tj2  ■  8nbo2 


-  Real 


lh'  (1 )' 

s 


(6.1) 


from  equation  (5.61),  It  is  therefore  possible  to  compute  the  despin 
moment  on  a  control  volume  of  dimensionless  length  2b  as  well  as  the  de¬ 
spin  moment  per  unit  length.  The  important  features  observed  from  Figure 


Pi 


»W  VW.-.'.N'y 


■f"  vSt ; 
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6  are,  firstly,  at  low  Reynolds  numbers  the  despin  moment  is  larger  in 
the  viscoelastic  cases  than  in  the  Newtonian  case;  secondly,  at  large 
Reynolds  numbers  the  effect  of  viscoelasticity  is  to  reduce  the  despin 
moment  by  comparison  with  the  Newtonian  case;  and  thirdly,  the  maximum 
despin  moment  occurs  at  a  lower  Reynolds  number  (Re  «  5-10)  for  the  vis¬ 
coelastic  liquids  than  for  the  Newtonian  liquid  (Re  =  15) .  These  suggest 
that  there  can  be  quite  significant  non-Newtonian  effects  in  the  precise 
low  Reynolds  number  regime  where  the  high  viscosity  instability  is  of 
practical  concern. 

It  is  interesting  to  compare  these  results  with  their  equivalents  for 
the  second-order  fluid,  for  which  S  is  given  by  (5.29).  In  Figures  7-10 
we  show  the  variation  of  f  and  g  with  radius  at  Reynolds  numbera  Re  -  15 
(Figures  7  and  8)  and  Re  -  50  (Figures  9  and  10).  In  addition  to  the 
Newtonian  curves  (Included  for  purposes  of  comparison)  we  show  the 
curves  for  values  of  Deborah  number  De  -  0.01  and  De  ~  0.1.  It  is  clear 
that  the  distortions  due  to  non-Newtonian  effects  are  very  different 
from  those  in  the  differential  model.  Figure  11  shows  the  quantity  M, 
defined  by  (6.1),  as  a  function  of  Reynolds  number  for  the  same  para¬ 
meter  values  De  -  0.01  and  De  -  0.1.  Again  the  differences  between  the 
differential  model  and  the  second-order  fluid  model  are  striking.  In  the 
latter  case,  as  shown  in  Figure  11,  although  the  moment  increases  due  to 
non-Newtonian  effects,  the  Reynolds  number  at  which  it  attains  its  peak 
is  greater  then  in  the  Newtonian  case. 

These  results  highlight  the  importance  of  a  proper  representation  of 
non-Newtonian  behavior  through  the  choice  of  a  constitutive  model  des- 


cribtng  the  fluid  motion  of  interest.  In  part  the  differences  between 
the  moments  shown  in  Figures  6  and  11  can  be  understood  from  the  follow¬ 
ing  argument.  The  value  of  the  moment  is  effectively  determined  by  the 
quantity  S  and  how  the  Bessel  functions  vary  with  S.  For  the  second- 
order  fluid,  from  (5.29), 


|  S  |  -  tt~— TTT  <  Re 
1  1  |l  +  lWe| 


(6.2) 


while  for  the  differential  model  (5.30) 


,  .  1  -  iWe 

S  “  Re  t: - - T  >  Re 

11  1  -  i<We 


(6.3) 


for  any  value  of  the  Weissenberg  number.  Thus  if  we  tentatively  take  the 
view  that  |s|  is  an  effective  Reynolds  number,  we  see  that  the  Reynolds 
number  is  reduced  (viscosity  is  increased;  shear  thickening)  for  the 
second-order  fluid,  while  the  Reynolds  number  is  Increased  (viscosity  is 
decreased;  shear  thinning)  for  the  differential  model.  The  general 
shapes  of  the  curves  in  Figures  6  and  11  are  consistent  with  these 
notions . 
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7.  FINITE  ELEMENT  METHOD 

The  problem  as  posed  in  Section  3  was  solved  numerically  by  a  finite 
element  method,  using  the  software  package  FIDAP  (Fluid  Dynamics  Analy¬ 
sis  Package).  FIDAP  is  a  commercial,  general  purpose  code  for  the  solu¬ 
tion  of  incompressible  fluid  flow  problems  governed  by  the  Navier-Stokes 
equations,  and  including  non- isothermal  effects.  The  code  applies  to 
both  transient  and  steady  flows,  and  can  handle  three-dimensional  prob¬ 
lems.  It  has  the  capability  of  addressing  flows  of  non-Newtonian  fluids. 

In  the  present  work  the  problem  was  posed  in  the  aeroballistic  reference 
frame  discussed  previously.  In  this  frame  the  problem  is  three- 
dimensional  and  a  steady  state  solution  is  required.  No  thermal  effects 
are  considered. 

The  essential  features  of  the  finite  element  method  as  incorporated  in 
FIDAP  are  as  follows.  The  domain  of  interest,  the  interior  of  the  cylin¬ 
der,  is  divided  into  a  number  of  geometrically  simple  elements,  thereby 
generating  the  finite  element  mesh  Since  the  code  FIDAP  works  in  carte¬ 
sian  coordinate  systems,  the  elements  are  straight-sided  even  though  a 
portion  of  the  cylinder  boundary  is  curved.  The  number  of  elements  used 
is  a  crucial  factor  in  determining  the  running  time  of  the  program,  and 
this  consideration  has  to  be  balanced  against  accuracy  requirements.  The 
finer  the  mesh,  the  greater  the  accuracy  and,  simultaneously,  the  cost; 
and  conversely,  coarsening  the  mesh  decreases  both  accuracy  and  cost. 
The  principal  criterion  in  deciding  on  the  mesh  is  that  the  associated 
approximations  should  be  sufficient  to  resolve  the  flow  field  every¬ 
where  ,  This  means,  for  example,  that  the  mesh  needs  In  he  quite  fine  in 
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regions  where  boundary  layers  occur  so  that  large  gradients  of  the  velo¬ 
city  field  can  be  appropriately  taken  into  account.  Thus,  flows  at  low 
Reynolds  numbers  can  generally  be  handled  with  a  fairly  coarse  mesh, 
while  high  Reynolds  number  flows  require  a  much  finer  mesh. 

In  the  present  project,  after  considerable  experimentation,  a  mesh  was 
designed  that  was  adequate  for  low  and  moderate  Reynolds  numbers,  over 
the  range  0  <  Re  <  2500  approximately.  At  the  lower  end  of  the  Reynolds 
number  range  a  coarser  mesh  would  no  doubt  have  sufficed,  but  it  was 
decided  to  retain  the  rather  fine  mesh  throughout  the  computations  in 
the  interests  of  obtaining  better  accuracy.  It  should  be  pointed  out 
that  the  Reynolds  number  referred  to  here  is  the  zero-shear-rate 
Reynolds  number;  in  the  non-Newtonian  situation  the  actual  Reynolds  num¬ 
ber  is  considerably  higher  due  to  shear  thinning. 

The  working  mesh  had  2240  elements,  which  were  8 -node  bricks,  and  2541 
nodes.  The  mesh  is  shown  in  Figure  12. 

In  the  finite  element  method  the  velocity  vector  is  approximated  on  each 
element  by  a  simple  polynomial  function.  In  the  present  case  the  veloci¬ 
ties  were  approximated  by  trilinear  interpolation  functions,  while  the 
pressure  was  represented  by  a  piecewise  constant  discontinuous  function. 
The  Galerxin  method  of  weighted  residuals  reduces  the  Navier-Stokes  and 
continuity  equations,  together  with  the  boundary  conditions,  to  a  large 
system  of  nonlinear  algebraic  equations  which  then  need  to  be  solved  by 
an  appropriate  technique.  The  solution  procedure  used  here  was  a  quasi- 
Newton  method,  which  is  an  iterative  method  of  Newton  type  involving  up¬ 
dating  of  the  iteration  matrix. 
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The  procedure  adopted  in  solving  the  problem  with  FIDAP  can  be  summar¬ 
ized  as  follows: 

1.  An  input  file  was  prepared  which  included, 

(a)  the  specification  of  the  mesh; 

(b)  the  specification  of  the  boundary  conditions,  namely  that 
the  container  is  in  rigid  body  rotation  in  the  aeroballistic 
frame,  by  a  call  to  a  subroutine  that  contains  the  boundary 
conditions  information; 

(c)  the  specification  of  the  solution  procedure,  chosen  from 
among  several  options  available  in  FIDAP; 

(d)  the  specification  of  the  non-Newtonian  properties  of  the 
liquid;  this  will  be  described  in  more  detail  in  the  next 
Section. 

2.  The  program  FIDAP  was  run  on  a  Hewlett  Packard  9000  computer,  it 
was  found  that  only  3-4  iterations  of  the  quasi-Newton  method  were 
required  to  achieve  convergence,  and  the  typical  run  time  on  this 
machine  was  about  12  hours  of  cpu  time.  The  first  run  was  at  a  low 
Reynolds  number,  and  subsequent  runs  at  higher  Reynolds  numbers 
were  performed  by  restarting  from  the  solution  obtained  in  the 
previous  run. 

3.  The  output  was  sent  to  a  post-processor  file  from  which  graphical 
information  concerning  the  flow  field  could  be  extracted. 
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4.  A  subroutine  was  used  to  compute  the  three  components  of  the 
moment  in  accordance  with  the  formula  (4.9). 


In  the  next  Section  we  detail  the  results  of  the  calculations. 
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8 .  COMPUTATIONS 

To  begin  with,  in  order  to  test  the  code  and  the  viability  of  the  mesh, 
we  performed  a  sequence  of  computations  for  a  Newtonian  liquid  over  a 
range  of  Reynolds  numbers.  The  aspect  ratio  of  the  cylinder  in  this  case 
was  b  -  4.368,  the  coning  angle  was  8  -  20°,  and  the  spin  ratio  was  r?  - 
1/6.  These  quantities  were  held  fixed  and  the  Reynolds  number  was  varied 
over  the  range  Re  -  10  to  Re  -  1000.  Figure  13  shows  the  variation  of 
the  roll  moment  Mz  with  Reynolds  number  over  this  range.  It  should  be 
noted  that  the  quantity  M  plotted  here  and  in  all  subsequent  graphs  is 
the  dimensionless  roll  moment  as  defined  by  equation  (4.9).  In  accord¬ 
ance  with  formula  (4.6)  the  conversion  to  dimensional  form  is  effected 
by  the  transformation 

M*  -  paia>ZMz  .  (8.1) 

It  can  be  seen  from  Figure  13  that  the  maximum  roll  moment  occurs  at  a 
Reynolds  number  approximately  equal  to  40.  This  is  larger  than  the  value 
predicted  by  Herbert  (1985)  but  the  latter,  of  course,  applied  to  the 
case  of  an  infinite  cylinder.  The  general  shape  of  the  curve,  however, 
is  consistent  with  that  obtained  from  the  infinite  cylinder  approxima¬ 
tion  as  shown,  for  example,  in  Figure  6. 

The  detailed  values  of  the  three  components  of  the  moment  are  shown  in 
Table  1.  The  Reynolds  numbers  listed  in  the  first  column  are  those  actu¬ 
ally  used  in  stepping  the  solutions  by  restarting  from  the  solution  of 
the  previous  run.  The  quantities  Mx,  Mz  ate  the  dimensionless  values 
defined  by  equation  (4.9),  while  the  quantity  is  the  pitching  moment 


due  to  deviation  from  solid  body  rotation,  namely 
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rt  -  My  -  rjoirb  (8.2) 

In  column  5  of  Table  1  we  have  written  the  ratio  Mz/Hx.  According  to  the 
calculation  presented  in  Section  4  this  ratio  should  exactly  equal 
tan  6  -  .3640  when  8  -  20°.  The  closeness  of  the  number  in  column  5  to 
this  value  is  one  good  indication  of  the  accuracy  of  the  computation. 

Typical  flow  fields  from  this  computation  are  shown  in  Figures  14-18. 
Figures  14 (a) , (b) , (c)  show  the  velocity  vector  in  the  transverse  plane 
defined  by  z  -  4.0  for  values  of  the  Reynolds  number  Re  -  20,  300,  1000 
respectively,  while  Figures  15(a) , (b) , (c)  show  the  equivalent  velocity 
vectors  at  the  same  respective  Reynolds  numbers  in  the  plane  z  -  -4.0. 
These  two  planes,  it  should  be  recalled,  are  quite  close  to  the  cylinder 
ends.  Figures  16 (a) , (b) , (c)  illustrate  the  contours  of  constant  axial 
vr  ocity  in  the  plane  z  -  4.0.  The  legends  in  the  Figures  show  the  reg¬ 
ions  of  upflow  and  downflow.  An  especially  interesting  feature  of  these 
diagrams  is  the  development  of  the  boundary  layers  at  the  wall  of  the 
cylinder;  the  boundary  layers  are  already  in  evidence  when  Re  -  300  and 
are  quite  strong  when  Re  -  1000.  Figures  17 (a) , (b) , (c)  show  the  ’’elocity 
vector  in  the  axial  cross-section  of  the  cylinder  lying  in  the  plane 
formed  by  the  vectors  w  and  0,  that  is,  the  plane  4>  -  0,  while  Figures 
18 (a) ,  (b) , (c)  depict  the  same  velocity  field  in  the  orthogonal  axial 
plane,  namely  <f>  -  n/2.  These  Figures  clearly  show  the  formation  of  two 
regions  of  strong  vortex- like  motion  as  the  Reynolds  number  increases. 

Next  we  performed  the  calculations  for  two  non-Newtonian  liquids  of 
interest  to  CRDC .  The  properties  of  these  liquids  were  presented  in  the 


form  of  experimental  data  in  which  the  variation  of  viscosity  with  shear 
rate  had  been  measured.  Some  data  on  the  variation  of  normal  stress  with 
shear  rate  was  also  available,  but  only  for  a  very  limited  range  of 
shear  rates.  A  constitutive  relation,  for  example  of  differential  type, 
was  not  known  for  either  of  these  liquids,  nor  was  it  possible  to  con¬ 
struct  one  on  the  basis  of  the  available  data.  The  measurements,  more¬ 
over,  were  understood  to  have  been  performed  in  simple  shearing  experi¬ 
ments,  in  which  the  nature  of  the  flow  was  very  different  from  that  in 
the  spinning  and  nutating  cylinder. 

In  the  absence  of  appropriate  data  for  the  liquids  under  the  flow  condi¬ 
tions  prevailing  in  the  present  problem,  it  was  decided  to  simulate  the 
non-Newtonian  effects  by  entering  a  shear- rate -dependent  viscosity  func¬ 
tion  extrapolated  from  the  data  provided.  This  data  is  shown  in  Table  2. 
One  of  the  liquids,  henceforth  referred  to  as  Liquid  1,  was  a  relatively 
low  viscosity  liquid,  with  a  zero-shear-rate  value  of  about  8  poise, 
wtiile  the  other,  Liquid  2,  was  a  high  viscosity  liquid  with  =  1300 
poise . 

The  calculations  were  performed  for  a  cylinder  of  aspect  ratio  b  -  4.50 
and  for  various  values  of  tj ,  P  and  Re,  This  aspect  ratio  applied  to  a 
test  fixture  having  radius  5.54  cm  and  length  49.84  cm.  Before  present¬ 
ing  the  results  of  the  computations,  we  note  that  the  Reynolds  number 
for  Liquid  1  corresponding  to  the  given  zero-shear-rate  viscosity  and  a 
spin  rate  of  4000  rpm  is  approximately  Re  -  1700,  and  increases  to 
Re  -  2500  for  spin  rate  of  6000  rpm.  For  Liquid  1,  therefore,  the  calcu¬ 
lations  are  close  to  the  limits  of  the  chosen  mesh.  For  Liquid  2,  the 
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corresponding  (Newtonian)  Reynolds  numbers  range  from  Re  «  10  for 
w  —  4000  r.p.ro.  to  Re  =  15  for  w  -  6000  r.p.m. 

Figures  19  and  20  relate  to  Liquid  1.  Figure  19  shows  the  variation  of 
dimensionless  roll  moment  H  with  Reynolds  number  in  the  case  fl  -  400 
rpm,  w  -  4000  rpm,  8  -  20°.  The  solid  curve  in  this  Figure  depicts  the 
behavior  for  the  given  non-Newtonian  Liquid  1,  while  the  broken  curve 
represents  the  behavior  of  the  "equivalent"  Newtonian  liquid,  that  is, 
the  Newtonian  liquid  having  viscosity  equal  to  the  zero-shear-rate 
viscosity  of  Liquid  1.  Figure  20  shows  a  repeat  of  these  calculations 
for  w  —  4500  rpm,  other  quantities  remaining  the  same.  We  see  that  at 
these  relatively  high  Reynolds  numbers  the  non-Newtonian  liquid  exper¬ 
iences  a  smaller  despin  moment  than  its  Newtonian  counterpart.  This 
result  is  consistent  with  the  predictions  of  the  analytic  approximation 
for  a  viscoelastic  fluid  as  shown  in  Figure  6. 

Figures  21  and  22  present  the  same  type  of  information  for  Liquid  2. 
Figure  21  shows  the  variation  of  with  Re  for  0  -  400  rpm,  u  -  4000 
rpm,  9  —  20°,  while  Figure  2  relates  to  the  spin  frequency  u>  -  4300  rpm. 
In  both  Figures  we  see  that  the  maximum  despin  moment  for  the  non- 
Newtonian  liquid  (solid  curves)  is  attained  at  very  low  Reynolds  num¬ 
bers,  much  lower  than  for  the  Newtonian  case.  Although  the  peak  values 
are  not  as  great  as  in  Figure  6,  the  general  picture  is  consistent  with 
that  result. 

Figure  23  also  relates  to  Liquid  2  and  shows  the  variation  of  ti  with 
spin  rate  for  some  different  values  of  coning  angle  and  coning  rate.  The 
curve  labelled  1  relates  to  the  case  0  —  500  rpm,  8  =  20°,  the  curve 
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labelled  2  relates  to  the  case  0  -  300  rpm,  8  -  20°;  and  the  curve  lab¬ 
elled  3  relates  to  the  case  0  -  500  rpm,  9  -  10°.  The  following  facts 
are  apparent  from  a  study  of  these  Figures:  the  roll  moment  decreases 
with  increasing  w  with  the  other  quantities  held  fixed;  at  fixed  w  and 
fixed  8  an  increase  in  coning  rate  from  300  rpm  to  500  rpm  results  in  an 
increase  in  Mz;  at  fixed  w  and  fixed  ft  an  increase  in  8  from  10°  to  20° 
results  in  an  increase  in  H  . 

X 

I  r  the  purposes  of  detailed  analysis,  the  data  shown  in  Figure  23  is 
presented  for  all  three  components  of  the  moment  in  Tables  3-5. 


9.  DISCUSSION 


In  reviewing  the  work  performed  in  this  project  we  note  several  salient 
features  that  merit  further  discussion. 

The  first  issue  is  the  relationship  between  the  results  obtained  by  the 
approximation  technique  of  Section  5  and  the  finite  element  computations 
of  Section  7.  In  discussing  this  question  we  leave  aside  temporarily 
whether  the  fluid  is  Newtonian  or  non-Newtonian,  since  the  question 
arises  equally  in  both  cases. 

More  specifically  the  point  is,  how  accurate  are  the  results  obtained  by 
the  methods  of  Section  5  as  presented  in  action  6?  The  procedure  used 
in  Section  5  involves  two  distinct  approximations  (and  consequently  two 
distinct  potential  sources  of  error):  it  is  assumed  that  the  cylinder  is 
infinitely  long,  and  it  is  assumed  that  the  coning  rate/spin  rate  ratio 
Is  so  small  as  to  allow  a  perturbation  in  powers  of  this  ratio,  »j ,  with 
only  the  leading  order  terms  retained.  It  should  be  emphasized  that 
these  two  assumptions  are  completely  independent,  and  a  solution  could 
formally  be  obtained  with  only  one  of  these  being  invoked.  The  infinite  - 
length  assumpt<or,  as  shown  in  Section  4,  directly  invalidates  any  calc¬ 
ulation  of  the  components  Mx  and  My  of  the  moment  because  of  the  loss  of 
rigid  end-wall  conditions,  but  leaves  the  possibility  of  a  sufficiently 
accurate  result  for  the  despin  moment  .  The  assumption  that  n  is  small 
on  the  other  hand,  introduces  an  error  in  the  determination  of  the  velo¬ 
cities  and  pressure,  and  hence  of  all  the  moments,  but  this  error  tends 
i.o  zero  as  17  -*  0. 


-49- 


The  finite  element  computation  gives  an  accurate  solution  provided  only 
that  the  physical  domain  is  sufficiently  well  resolved,  that  is,  that 
the  mesh  is  fine  enough.  The  mesh  used  in  our  calculations  was  tested 
extensively  and  found  to  be  adequate  for  the  range  of  Reynolds  numbers 
under  consideration.  Therefore  we  regard  our  computational  results  as 
highly  accurate,  especially  at  lower  Reynolds  numbers. 

As  one  indication  of  the  validity  of  the  infinite-cylinder,  small  rj  per¬ 
turbation  solution,  we  chow  in  Table  6  values  for  the  despin  moment 
calculated  by  the  two  methods.  The  results  refer  to  a  Newtonian  liquid, 
with  b  -  4.368,  6  -  20°  and  rj  -  1/6.  The  perturbation  results  are  taken 
from  Figure  6,  converted  using  the  formula  (6.1),  while  the  computa¬ 
tional  results  are  those  shown  in  Figure  13.  We  see  from  this  Table  that 
the  values  differ  by  between  5%  and  25%,  depending  on  Reynolds  number, 
and  that  the  maximum  of  the  despin  moment  occurs  at  about  Re  -  20  in  one 
case  and  at  about  Re  -  40  in  the  other. 

The  second  major  issue  is  that  of  the  appropriate  representation  of  the 
non-Newtonian  behavior  of  a  liquid.  There  are  two  aspects  to  this: 
selection  of  a  suitable  constitutive  relation  and  determination  of  the 
parameters  to  be  inserted  in  this  relation.  It  is  crucial,  moreover, 
that  the  relation  and  the  parameters  be  appropriate  for  the  type  of  flow 
under  consideration.  In  our  view,  a  constitutive  relation,  with  para¬ 
meters,  that  has  been  determined  and  tested  f<-  iifferent  type  of  flow 
may  well  be  completely  irrelevant  for  some  ot  -pe  of  flow. 

This  situation  poses  a  major  dilemma  for  the  present  project.  On  a  gen¬ 
eral  level  we  are  unaware  of  any  constitutive  relations  that  have  been 
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verified  (as  against  postulated)  for  fully  three-dimensional  flows  of 
the  type  under  study  here.  For  this  reason  we  were  not  in  a  position  to 
hypothesize  a  constitutive  law  that  could  be  subject  to  experimental 
comparison.  On  a  specific  level,  for  the  two  liquids  of  direct  interest 
the  information  available  regarding  their  non-Newtonian  behavior  was 
obtained  from  simple -shear  experiments,  and  may  or  may  not  be  relevant 
to  the  flow  in  a  spinning  and  coning  cylinder. 

As  a  consequence  of  this  we  cannot  be  certain  whether  the  results 
obtained  in  our  computations  have  quantitative  validity.  On  the  other 
hand  we  feel  that  they  do  predict  correct  behavior  qualitatively;  this 
is  because  there  is  general  qualitative  agreement  between  the  computa¬ 
tions  based  on  the  experimental  data  and  the  approximate  analysis  in 
Section  5  based  on  a  number  of  different  theoretical  models.  The  most 
important  feature  that  emerges  uniformly  is  that  the  maximum  despin 
moment  occurs  at  a  significantly  lower  Reynolds  number  in  a  non- 
Newtonian  fluid  than  in  a  Newtonian  fluid. 

In  order  to  obtain  computational  results  that  could  be  quantitatively 
reliable  it  would  be  necessary  to  determine  a  constitutive  law  and  para¬ 
meters  valid  for  the  flow  in  question.  This  would  not  be  easy  to  do  with 
any  degree  of  certainty.  The  analysis  of  Section  5  can  be  helpful;  it 
shows  that  beginning  with  various  different  models  one  arrives  for  small 
enough  rj  at  a  relation  that  involves  a  complex  viscosity.  The  latter  is 
a  quantity  that  can  be  measured  from  stress  relaxation  or  torsional 
rheogoniometer  experiments.  At  the  same  time  it  must  be  said  that  a  com¬ 
plex  viscosity  law,  which  is  a  linear  law,  may  be  Insufficient.  The  non- 


linear  law  of  shear  chinning  used  in  the  computations  led  to  significant 
effects,  which  would  not  have  been  the  case  if  the  response  of  the 
liquid  were  primarily  linear-elastic.  This  nay  be  attributed  to  the  fol¬ 
lowing  fact.  A  typical  dilute  polymer  solution  has  a  relaxation  time  of 
about  0.05  -  0.1  sec.  For  a  large  spin  rate  of  4000-6000  rpm,  the 
Weissenberg  number  as  defined  by  equation  (2.6)  will  lie  in  the  range 
20-b0,  which  is  very  large  in  terms  of  standard  viscoelastic  models.  In 
fact,  the  perturbation  procedure  of  Section  5  breaks  down  at  such  large 
Weissenberg  numbers.  These  observations  suggest  an  appropriate  empirical 
law  should  incorporate  both  the  elastic  effects  giving  rise  to  complex 
viscosity  and  nonlinear  effects  due  to  the  large  Weissenberg  number. 

In  summary  therefore,  we  believe  that  the  finite  element  method  as 
implemented  in  the  code  FIDAP  provides  accurate  solutions  to  the  problem 
of  non-Newtonian  fluid  flow  in  a  coning  and  spinning  cylinder;  reliable 
quantitative  results  await  the  empirical  determination  of  a  suitable 


constitutive  law. 
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Re 

M 

X 

H 

7 

M 

t 

M  /M 

Z/  X 

10 

.05953 

.02807 

.02202 

.3699 

20 

.07791 

.06560 

.02860 

.3671 

30 

.08056 

.09023 

.02946 

.3657 

40 

.07982 

.10545 

.02918 

.3656 

60 

.07654 

.12556 

.02805 

.3665 

80 

.07303 

.13854 

.02681 

.3671 

100 

.06975 

.14793 

.02562 

.3673 

120 

.06669 

.15532 

.02450 

.3674 

150 

.06288 

.16327 

.02311 

.3675 

175 

.05984 

.16885 

.02200 

.3676 

200 

.05729 

.17310 

.02108 

.3680 

300 

.04921 

.18448 

.01817 

.3692 

400 

.04361 

.19069 

.01619 

.3712 

500 

.03945 

.19437 

.01472 

.3731 

750 

.03244 

.19805 

.01229 

.3789 

1000 

.02804 

.19768 

.01080 

.3852 

Table  1:  Moments  for  a  Newtonian  liquid  in  a  cylinder  of  aspect 
ratio  4.368,  coning  angle  20°,  coning/spin  ratio  1/6 


Liquid  l 


Liquid  2 


Shear  rate 
(sec*1) 

Viscosity 

(poise) 

Sheer  rate 
(s-c'1) 

Viscosity 
(poise ) 

.855 

8.011 

.017 

1307.349 

1.707 

7.937 

.034 

1216.849 

3.407 

7.488 

.067 

1125.913 

6.798 

7.271 

.135 

1022.781 

13.564 

6.818 

CM 

942.722 

27.064 

5.892 

.54 

797.31 

54. 

4.872 

1.077 

651.202 

107.744 

4.07 

2.149 

504.397 

214.977 

3.152 

4.289 

371.76 

428.937 

2.49 

8.558 

257.124 

679.819 

2.087 

17.076 

161.551 

Table  2:  Experimental  data  for  variation  of  viscosity 
with  shear  rate  for  two  liquids 


Table  3:  LIQUID  2  -  Coning  rate  0  -  500  rpm,  Coning  angle  0  =  20 


rpm 

rj 

Re 

M 

X 

"y 

H 

z 

1000 

.3 

1.7 

.23166 

.113717 

.083304 

2000 

.15 

3.4 

.067447 

.071241 

.024475 

3000 

.1 

5.1 

.028997 

.040781 

.010542 

4000 

.075 

6.8 

.015686 

.026023 

.005706 

5000 

.06 

8.5 

.009542 

.017878 

.003479 

6000 

.05 

10.2 

.0064163 

.0131311 

.0023306 

7000 

.04286 

11.9 

.004534 

.0099828 

.0016464 

8000 

.0375 

13.6 

.0033457 

.0078472 

.0012149 

9000 

.0337 

15.3 

.0025481 

.00631841 

.00092651 

10000 

.  -3 

17.0 

.002003 

.005196 

.000728 

Table  4:  LIQUID  2  -  Coning  rate  fi  -  300  rpm,  Coning  angle  6  -  20° 


u> 

rpm 

Re 

H 

X 

M 

y 

M 

t 

1000 

.5 

1.7 

.34027 

.160873 

.059135 

2000 

.25 

3.4 

.10254 

.103172 

.017882 

3000 

.1667 

5.1 

.043883 

.0600640 

.0077028 

4000 

.125 

6.8 

.023523 

.0383475 

.0041432 

5000 

.1 

8.5 

.014343 

.0266239 

.0025288 

0000 

.0833 

10.2 

.0095728 

.0195625 

.0016808 

7000 

.07143 

11.9 

.0067004 

.0149356 

.0011818 

6000 

.0625 

13.6 

.0049337 

.0117658 

.000869 

9000 

1 

.05556 

15.3 

.0037631 

.009499 

.000663 

10000 

.05 

17.0 

.0029458 

.007813 

.000519 

Table  5:  LIQUID  2  -  Coning  rate  0  -  500  rpra,  Coning  tingle  t)  -  10 


Re 

M 

z 

finite  cylinder 
computation 

r 

M 

Z 

infinite  cylinder 
approximate  solution 

10 

.02202 

.02656 

20 

.02860 

.03043 

30 

.02946 

.02882 

40 

.02918 

.02690 

60 

.02805 

.02386 

80 

.02681 

.02268 

100 

.02562 

.01998 

120 

.02450 

.01865 

ISO 

.02311 

.01709 

200 

.02108 

.01520 

Table  6:  Comparison  of  despin  moments  obtained  from  finite  element 
computation  and  infinite-cylinder  approximate  solution 
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Figure  6:  Variati 
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Figure  7:  Variation  of  velocity  function  f  with  radius  r 
for  second-order  fluid,  Re  -  15 
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Figure  8:  Variation  of  velocity  function  g  with  radius  r 
for  second-order  fluid,  Re  -  18 

1.  Newtonian 

2.  De  -  .01 

3.  De  -  .1 
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Figure  14c 
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SPINNING  S  NUTATING  CYLINDER  :  Re  -  300 
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Variation  of  despin  moment  Mz  with  Reynolds  number 
for  non-Newtonian  Liquid  2  (solid  curve)  and  for  an 
equivalent  Newtonian  liquid  (broken  curve)  with 
Q  -  400  rpm,  w  -  4000  rpm,  5-20° 
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Figure  22:  Variation  of  despin  moment  Mz  with  Reynolds  number 
for  non-Newtonian  Liquid  2  (solid  curve)  and  for  an 
equivalent  Newtonian  liquid  (broken  curve)  with 
-  Q  =  400  rpm,  w  =  4800  rpm,  9  =  20° 


